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Abstract. We show that the well-known discrepancy, known for about two decades, between the radial dependence of the 
Galactic cosmic ray nucleon distribution, as inferred most recently from EGRET observations of diffuse 7-rays above 100 MeV, 
and of the most likely cosmic ray source distribution (supernova remnants, superbubbles, pulsars) can be explained purely by 
propagation effects. Contrary to previous claims, we demonstrate that this is possible, if the dynamical coupling between the 
escaping cosmic rays and the thermal plasma is taken into account, and thus a self-consistent calculation of a Galactic Wind is 
carried out. Given a dependence of the cosmic ray source distribution on Galactocentric radius r, our numerical wind solutions 
show that the cosmic ray outflow velocity, V(r, z) = uq + Vao, also depends both on r, as well as on vertical distance z, 
with uo and Vao denoting the thermal gas and the Alfven velocities, respectively, at a reference level zc- The latter is by 
definition the transition boundary from diffusion to advection dominated cosmic ray transport and is therefore also a function 
of r. In fact, the cosmic ray escape time averaged over particle energies decreases with increasing cosmic ray source strength. 
Thus an increase in cosmic ray source strength is counteracted by a reduced average cosmic ray residence time in the gas 
disk. This means that pronounced peaks in the radial distribution of the source strength result in mild radial 7-ray gradients at 
GeV energies, as it has been observed. The effect might be enhanced by anisotropic diffusion, assuming different radial and 
vertical diffusion coefficients. In order to better understand the mechanism described, we have calculated analytic solutions of 
the stationary diffusion-advection equation, including anisotropic diffusion in an axisymmetric geometry, for a given cosmic 
ray source distribution and a realistic outflow velocity field V{r, z), as inferred from the self-consistent numerical Galactic 
Wind simulations performed simultaneously. At TeV energies the 7-rays from the sources themselves are expected to dominate 
the observed "diffuse" flux from the disk. Its observation should therefore allow an empirical test of the theory presented. 
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1. Introduction 

Our information on the spatial distribution of cosmic rays 
(CRs) in the Galaxy stems largely from measurements of non- 
thermal emission, generated by the energetic charged parti- 
cles interacting with matter and electromagnetic fields. For 7- 
ray energies above 100 MeV, the main production process is 
probably via 7r"-decay, resulting from nuclear collisions be- 
tween high energy particles and interstellar matter. Past and 
recent observations in the GeV range have shown a roughly 
uniform distribution of diffuse 7-ray emissivity in the Galactic 
plane, exhibiting only a shallow radial gradient (in a cylindri- 
cal Coordinate system). Hence, if 7-rays were to map the spa- 
tial CR distribution, we would expect it to be uniform as well. 
However, associating CR production regions with star forma- 
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tion regions, all possible Galactic CR source distributions are 
strongly peaked towards a Galactocentric distance at which a 
ring of molecular gas resides. It is commonly believed that the 
bulk of the CR nucleons below about 10^^ cV is produced in 
supernova remnants (SNRs) - the majority being core collapse 
SNRs - most likely by diffusive shock acceleration (Krymsky 
1977; Axford et al. 1977; Bell 1978a,b; Blandford & Ostriker 
1978). It has been argued that up to some tenths of the hydro- 
dynamic explosion energy might be converted into CR energy 
(e.g., Berezhko et al. 1994). Since high-mass star formation 
mostly occurs in a spatially nonuniform manner, i.e. in OB as- 
sociations predominantly located in the spiral arms of late type 
galaxies, we are confronted with the problem of reproducing 
a mild radial gradient in the diffuse Galactic 7-ray emission 
as it has been observed for the first time by the COS-B satellite 
(Strong et al. 1988) and, more recently, with higher angular and 
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energy resolution, higher sensitivity and lower background by 
the EGRET instrument of the CGRO satellite (Strong & Mattox 
1996, Digel et al. 1996). If the SNRs are the sources of the CR 
nucleon component and if the source distribution is inhomo- 
geneous, this discrepancy must arise during the propagation of 
CRs from their sources through the interstellar medium. Unlike 
e.g. the interpretation of radio synchrotron emission, gener- 
ated by relativistic electrons, 7-ray data open the possibiUty 
of studying the nucleonic component of the CRs, in which al- 
most all of the energy is stored (see e.g., Dogiel & Schonfelder 
1997). The distribution of the 7-ray emissivity in the Galactic 
disk therefore bears important information on the origin of CRs 
and on the conditions of CR propagation in the Galaxy. 

2. Gradient of the 7-Ray Emissivity in tlie 
Galactic Disk 

The first data on the radial distribution (i.e. gradient) of the 7- 
ray emissivity in the disk were obtained with the SAS-II satel- 
Ute (energy range 30 - 200 MeV). The data showed rather 
strong variations of the emissivity along the Galactic plane 
which dropped rapidly with radius (see e.g., Stecker & Jones 
1977), in rough agreement with the distribution of candidate 
CR sources in the Galaxy such as SNRs or pulsars. A notice- 
able discrepancy however emerged from the COS-B data (en- 
ergy range 70 - 5000 MeV), in which the emissivity gradient 
was found to be rather small, when compared to the SNR dis- 
tribution, especially at high energies (E^ >300 MeV), with 
a maximum variation by a factor of only 2 (see Strong et al. 
1988). An energy dependence of the gradient was mentioned 
by almost all groups analyzing the COS-B data and was usu- 
ally interpreted as a steeper gradient for CR electrons, produc- 
ing the soft part of the 7-ray spectrum, compared to nuclei. 
Detailed analyses of different models of CR propagation based 
on the gradient value were usually performed with the COS-B 
data because of their much better statistics compared with the 
SAS-II data. 

The recent measurements with the EGRET instrument on- 
board CGRO, obtained by different methods at energies 100 
- 10000 MeV, showed that the emissivity drops at the edge of 
the disk. Digel et al. (1996), performed a study of the outer part 
of the Galaxy towards molecular complexes, the Cepheus and 
Polaris Flares in the local arm and the larger molecular com- 
plex in the Perseus arm. Since the total masses of these com- 
plexes are known, it was possible to infer the CR density from 
the measured 7-ray fluxes in the direction of these objects. It 
was discovered that the apparent emissivity decreases by a fac- 
tor 1.7, which is somewhat smaller than that for COS-B. The 
analysis of Strong & Mattox (1996), based on a model of the 
average gas density distribution in the disk, shows a smaller in- 
tensity gradient which does not differ significantly from that of 
COS-B in the Galactic disk. 

One of the important conclusions which follows from all 
these data, and which we wiU discuss below, is that the distri- 
bution of the 7-emissivity in the disk in the GeV range is rather 
uniform compared with the most probable distribution of CR 
sources. 



3. Basic Ideas for a New Approach of CR 
Propagation in the Galaxy 

A natural explanation of a uniform CR distribution would be 

effective radial mixing due to the diffusion of CRs produced 
in different parts of the Galactic disk. It is then straightfor- 
ward to infer the mixing volume, which usually includes the 
Galactic disk plus a large Galactic halo; the details of such a 
model were discussed in the book of Ginzburg and Syrovatskii 
(1964). This 3-dimensional cyUndrical model, which includes 
diffusive-advective transport, a free escape of cosmic rays from 
the halo boundary into intergalactic space, and the observed su- 
pernova shells as sources of cosmic rays in the Galactic disk, 
explained very well characteristics of the Galactic radio and 
7-ray emission as well as the data on cosmic ray spectra and 
their chemical composition including stable and radioactive 
secondary nuclei, intensities of positrons and antiprotons etc. 
(a summary of this analysis can be found in Berezinskii et al. 
1990). 

In general the propagation equation for CRs is described by 

d f dE \ 

- V,(A,Vjn - U,n) + — f — nj = Qir,E) , (1) 

where Dij is the diffusion tensor and Ui is the advection ve- 
locity, which in general are functions of the coordinates, and 
dE/dt is the term which describes particle energy losses in- 
cluding adiabatic losses if the velocities vary spatially. The 
right hand side of the equation describes the spatial and en- 
ergy dependence of the source density, derived from radio ob- 
servations of the SNR distribution and the injection spectrum 
of CRs. 

The main conclusion was that almost all observations can 
be reasonably explained if the halo extension is about several 
kpc, the injection spectrum of electrons and protons Q{E) is a 
power-law (Q{E) oc E''''" with 70 ~ 2.2 - 2.4), and the radi- 
ally uniform velocity of advection does not exceed the value of 
20 km/s. Compton scattering of relativistic electrons was found 
to play a significant role in the halo 7-ray emission. However, 
this model failed to explain a smooth emissivity distribution of 
7-rays in the Galactic disk (see for details Dogiel & Uryson 
1988). Even in the case of a very extended halo with a radius 
larger than 10 kpc the derived emissivity gradient (calculated 
for the observed SN distribution) was larger than observed. 
Only for a hypothetical uniform distribution of the sources in 
the disk the calculations can reproduce the data. 

This model was further developed by Bloemen et al. (1991, 
1993). Extensive investigations of the 3-dimensional diffusion- 
advection transport equation for nucleons, low-energy elec- 
trons and the 7-rays showed that even in the most favorable 
case of an extended halo with a vertical height as large as 
20 kpc, their model, although reproducing the COS B data 
marginally, is not able to remove the signature of the observa- 
tionally inferred SNR source distribution, i.e. the distribution 
of the calculated emissivity was still steeper than permitted by 
the data. The vertical gradient of the advection velocity derived 
from the fluxes of stable and radioactive nuclei near Earth had 
to be smaUer than 15 km/s/kpc. The other unexpected conclu- 
sion was that the halo extension obtained from the nuclear data 
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was significantly less than estimated from the 7-ray data (see 
also Webber et al. 1992). 

Recently, Strong and collaborators in a series of papers 
(see Strong et al. 2000, and references therein) have devel- 
oped a numerical method for this model and made a new at- 
tempt to analyze the 7-ray emission and the CR data, based on 
the latest data from COMPTEL, EGRET and OSSE. They lim- 
ited the analysis to radially uniform diffusion-advection and to 
diffusion-reacceleration transport models, and concluded that 
no such diffusion-advection model can adequately describe the 
data, in particular the B/C ratio and the energy dependence. 
These authors claimed, however, that by including reaccelera- 
tion one can account for all the observational data. A peculiar 
consequence of their analysis was that they did not use the ob- 
served supernova source distribution (see e.g., Kodeira 1974, 
Leahy & Xinji 1989, Case & Bhattacharya 1996, 1998) as the 
input Q for Eq. (1), but rather derived it from the 7-ray data 
to reproduce the observed spatial variations of the emissivity in 
the disk. The result was a source distribution that is flat in the 
radial direction. 

Thus we see that the conventional model of unifrom 
diffusion-advection has serious problems in spite of its evident 
achievements. One solution is to assume that some of the obser- 
vational data are not significant like the SN distribution derived 
from the radio observations (Strong et al. 2000). The alterna- 
tive is to conclude that it is time to abandon the standard model, 
which is what we do in this paper We shall demonstrate, that 
strong radial source gradients will be removed by a strong ad- 
vection velocity in the halo (due to a Galactic wind driven by 
the CRs themselves, see below) that varies with radius R and 
height z. In addition anisotropic diffusion with different diffu- 
sion coefficients and D± in the disk and the halo, respec- 
tively, might also play a role. It should be emphasized that a 
radially varying advection velocity occurs naturally in spiral 
galaxies, even for a uniform source distribution, because the 
gravitational potential increases towards the centre, thus induc- 
ing stronger velocity gradients in this direction (see Fig. ||), as 
has been shown by Breitschwerdt et al. (1991). 

The existence of strong advective CR transport in the 
Galactic halo has been shown on dynamical grounds in a num- 
ber of papers in the past (Ipavich 1975; Breitschwerdt et al. 
1987, 1991, 1993;Fichtneretal. 1991; ZirakashviU et al. 1996; 
Ptuskin et al. 1997). The key element of halo transport theory is 
that CRs, which by observations are known to escape from the 
Galaxy, resonantly generate waves by the so-called streaming 
instability (Kulsrud & Pearce 1969) leading to strong scatter- 
ing of CRs. Therefore, even in the case of strong non-linear 
wave damping, advection is at least as important a CR trans- 
port mechanism out to large distances in the halo as diffusion, 
provided that the level of MHD-turbulence is high enough for 
coupling between CRs and MHD waves (Dogiel et al. 1994, 
Ptuskin et al. 1997). There is also growing indirect observa- 
tional evidence of outflows from the interpretation of soft X- 
ray data of galactic halos in edge-on galaxies like NGC 463 1 
(Wang et al. 1995), and also of the soft X-ray background in 
our own Galaxy (Breitschwerdt & Schmutzler 1994, 1999). 
Furthermore, the near constancy of the spectral index of non- 
thermal radio continuum emission over large distances along 



the minor axis in the halo of edge-on galaxies is most natu- 
rally explained by an advective transport velocity of relativistic 
electrons (along with the nucleons) that is ever increasing with 
distance from the Galactic plane (Breitschwerdt 1994). 

In the disk of spiral galaxies, the regular magnetic field is 
following roughly the spiral arms and is therefore mostly par- 
allel to the disk, with noticeable deviations in some regions 
where outflow is expected. Here, also a regular vertical compo- 
nent seems to be present (eg. Hummel et al. 1988), which has 
been detected in a number of galaxies like NGC 4631, NGC 
5775 (Tullman et al. 2000), and NGC 4217. Multi-wavelength 
observations of the galaxy NGC 253 show a local correlation 
between non-thermal radio continuum. Ha and X-ray emission 
near the disk-halo interface in off-nuclear regions (Dettmar 
1992; M. Ehle, private communication). This also spatially co- 
incides with enhanced star formation activity in the disk as can 
be seen from FIR data. 

Since the disk is not fully ionized in contrast to the halo 
and since waves are efficiently dissipated there by ion-neutral 
damping, the most important contribution to the random field 
in the disk is by turbulent mass motions, induced by supernova 
explosions and other stellar mass loss activity. Thus the wave 
spectrum will be very different from the one in the halo, where 
self-excited waves subject to nonlinear wave damping (Dogiel 
et al. 1994, Ptuskin et al. 1997), satisfying the gyro-resonance 
condition, dominate. Consequently, the averaged diffusion co- 
efficient will be different in the Galactic disk and the halo (see 
Sect. ^). Therefore, we believe that anisotropic diffusion, to- 
gether with radially varying advection, is the most general and 
most probable mode of CR transport to occur. 

We investigate the CR transport processes of diffusion and 
advection and discuss the possibility of flattening radial CR 
source gradients of a given SNR distribution (despite observa- 
tional uncertainties) by particle propagation. In our view it is 
using the "natural" boundary condition when calculating the 
response of CR transport to a given CR source distribution 
(whatever its observational limitations may be at the present 
time). This is in contrast to adjusting the source distribution a 
posteriori. To implement the transport processes properly we 
shall allow for radially varying advection, anisotropic diffu- 
sion, (different values of the diffusion coefficient parallel and 
perpendicular to the disk) and the appropriate boundary condi- 
tions to employ a cosmic ray distribution and then to calculate 
from these the resulting spatial distribution of CRs. 

In Sect. ^ we discuss the question of anisotropic CR diffu- 
sion and in Sect. ^ we introduce a simple model that describes 
the observed CR source distribution. In the following Sect. ^ 
we show heuristically (and in Sect. || also analytically) how a 
radially varying CR source distribution induces variations in 
the CR energy density, which in turn leads to a radial varia- 
tion of the diffusion-advection boundary, zq, and the outflow 
velocity, respectively, and thus to a tendency to flatten the ra- 
dial CR distribution and hence the 7-ray emissivity gradient. 
In Sect. ^ we demonstrate by numerical calculations how such 
a source distribution will naturally generate a radially vary- 
ing outflow velocity. In Sect. || we discuss in detail models of 
different complexity for advective-diffusive transport with ra- 
dially varying outflow velocity, and show analytically how in 
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each case for a given CR source distribution its radial signa- 
ture on the energetic particle and 7-ray distribution is reduced. 
The most advanced of these models include a dependence of 
galactic wind velocity on the CR source strength. In addition a 
full analytic solution of the 3-dimensional diffusion-advection 
equation in axisymmetry for a given realistic velocity depen- 
dence u{r, z), parallel and perpendicular to the Galactic disk 
similar to the one derived in Sect, ^is calculated. In Sect. ^ the 
source contributions to the "diffuse" 7-ray background at TeV 
energies are taken into account and shown to be highly signifi- 
cant. This should allow an empirical test of our theoretical pic- 
ture. In Sect. |l^ we discuss and summarize our results. A num- 
ber of detailed calculations can be found in the Appendices ^ 
Handg 

4. Anisotropic Diffusion Coefficient 

In most models of diffusive CR propagation, the diffusion ten- 
sor is approximated by a scalar quantity D, representing spa- 
tially uniform transport, D = DqE"^, where a = 0.5 — 0.6 
describes the energy dependence. However, there are good rea- 
sons why the diffusion coefficient may be anisotropic. 

The propagation of CRs in the interstellar medium is 
mainly determined by their interaction with electrical and mag- 
netic fields. CRs interact strongly with fluctuations of the mag- 
netic field (MHD-waves) and are scattered by them in pitch 
angle. In the simplest case the total magnetic field consists of 
two components and can be written as 



B = Bq + 5B, 



(2) 



where Bq is the average large scale magnetic field and 5B is 
the field of small scale fluctuations. 

Effective scattering of particles by these fluctuations occurs 
when the interaction is resonant, i.e. when the scale of the fluc- 
tuations in the magnetic field Bq is of the order of the par- 
ticle gyroradius. This leads to a stochastic motion of particles 
through space; the associated diffusion coefficient D\\ along the 
magnetic field Bq is of the order of 



(3) 



where v is the particle velocity. The value L|| ~ v/v is the 
transport length of the particle along the magnetic field lines, 
and V is the rate of resonant particle scattering. The main as- 
sumption of this description is that the fluctuation amplitude is 
much less than the average magnetic field 



SB<^Bo 



(4) 



The displacement of a charged particle perpendicular to the 
regular magnetic field Bq due to scattering is therefore small: 



Du 



(5) 



I.e., 

^11 » , (6) 
where rg is the particle gyroradius, and Vg ^ Ln. 



The situation is more complex if there exists also a large 
scale random magnetic field AB whose scale is much larger 
than the particle gyroradius: 



B = Bo + AB + 5B . 



(7) 



In this case the diffusion transverse to the magnetic field Bq 
occurs due to the divergence of the particle trajectories, which 
is much faster than the resonant diffusion from Eq. || 

The procedure to derive the transport equation for CRs in 
this case was described , e.g. by Toptygin (1985), who showed 
that the maximum value of D± is: 



(8) 



Thus, diffusion in the Galaxy is quasi-isotropic only in the case 
AB - Bo. 

A more precise analysis (see Berezinsky et al. 1990, chap. 
9) shows, however, that in the interstellar medium the correla- 
tion between the components of the diffusion tensor leads to an 
effective perpendicular diffusion coefficient: 



Df - 0.1 • Dn 



(9) 



Observations show that in the Galactic disk the large scale 
magnetic field is typically parallel to the disk, whereas in the 
halo also a significant perpendicular component may exist, as 
it has been observed in the case of the edge-on Sc galaxy NGC 
4631. Therefore one can expect ^ Ki in the disk, and 
Dz ^ Ky in the haloQ with and being the components of 
the diffusion tensor perpendicular and parallel to the Galactic 
plane. Thus, in general, CRs spend a considerable amount of 
their transport time in the halo, subject to anisotropic diffusion. 

5. A simple modei for tlie Gaiactic CR source 
distribution 

The inference on the spatial distribution of CR sources from 
direct observations is plagued by a number of problems. From 
energy requirements for the bulk of CRs below 10^^ eV, it is 
known that the only non-hypothetical Galactic candidates for 
the sources are SNRs and pulsars, with global energetic re- 
quirements favouring the former SNRs can be best studied 
in the radio continuum and in soft X-rays, but as low surface 
brightness objects larger and older remnants are systematically 
missed. Since samples are usually flux limited, the more distant 
objects will be lost as well. Pulsars on the other hand should be 
present in the Galaxy in large numbers. However they are only 
detectable if their narrow beams happen to cross the line of 
sight of the observer or, in X-rays, as isolated old neutron stars; 
so far only four candidates are known from the ROSAT All- 
Sky Survey (Neuhauser & Triimper 1999). Therefore selection 
effects will bias samples heavily in both cases. 

Based on the observations of SNRs by Kodeira (1974) and 
pulsars by Seiradakis (1976), Stecker & Jones (1977) have 
given a simple radial distribution of the form 



Qi{x) = qox'" exp{-bx) , 



(10) 



' Note that in the disk, Hr, and in the halo, D^, are both parallel to 
the large scale regular magnetic field direction. 
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where, Qi{x) is a SNR surface density in x, with x = t/Rq, 
and a and b are matched to the observations. Here we adopt a = 
1.2 and b = 3.22. Since both SNRs and pulsars are associated 
with Population I objects, it is reassuring that a and b are very 
similar for both classes of objects. 

To determine the proportionality constant qa for the dis- 
tribution given in Eq. ([lO|), we use the results obtained by 
Leahy & Xinji (1989) who used the catalogue of Li (1985) 
derived from radio observations. Leahy & Xinji (1989) con- 
sidered shell-type SNRs and applied empirical correction fac- 
tors due to the incompleteness of the flux limited sample. The 
spatial distribution thus obtained shows a peak at 4-6 kpc from 
the Galactic center (see Fig. [l|), assuming an overall rotational 
symmetry. A more systematic study has recently been under- 
taken by Case & Bhattacharya (1998), who have made a care- 
ful analysis of an enhanced Galactic SNR sample, using an im- 
proved E — relation. These authors find a peak in the dis- 
tribution at around 5 kpc^ and a scale length of ^ 7 kpc, thus 
confirming the gross features of the older work. 

The number of SNRs in an annulus of width dx is given by 



dN = Qi{x)2Trxdx 

= 2Trqox"'^^ exp{—bx)dx , 



(11) 



and thus the total number of SNRs in the disk of radius 
18 kpc is 



Dragicevich et al. (1999) have analyzed also the radial dis- 
tribution of SNe in a sample of 218 external galaxies of dif- 
ferent Hubble types, corrected for inclination angle. The data 
were sorted into radial bins and the numbers converted into 
SNe surface densities. They showed that the radial SN surface 
distribution can be well fitted by an exponential radial decrease 
of the form 



S{R) = CTc e-^/^" 



(15) 



with (Tc = 350. Taking a sub-sample of 36 Sb-Sbc galaxies, 
they were able to derive a value of Rsc — 5.4 kpc for the ex- 
ponential scale length of the Milky Way as a representative Sb- 
Sbc galaxy. The number of SNe derived is consistent with his- 
torical records of Galactic SNe within 4 kpc from the Sun. In 
Fig. |l] it is shown that both the exponential decrease and, in- 
terestingly, but rather fortuitously, also the absolute numbers in 
the SN surface density agree fairly well for i? > 5 kpc with the 
earlier derived relation ( p^ for the SNR surface density. 

Therefore the value calculated for the proportionality con- 
stant, qq, provides a reliable estimate for the number of Galactic 
SNRs in a circular ring at a distance R. 

Radial SN Distribution in galaxies 



iV = 27r / qqx"-'^^ exp{-bx)dx 
Jo 



2nqo / X exp{—bx)dx 
Jo 
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for Xa — 1.8, awl and b = 3.22. Numerical integration of 
Eq. 1^ yields TV = 0.33(70. 

On the other hand, the total number of SNRs that should 
be observable at present is their production rate times their 
average life time tsn before they merge with the hot intercloud 
medium and lose their individual appearance, or observation- 
ally escape the detection limit. The rate at which SNe occur in 
the disk is 2.2 ± 1.3 per century for randomly exploding stars 
(van den Bergh 1990, 1991) and 0.45 per century for explo- 
sions occurring in OB associations (Evans et al. 1989), giving 
a total SN rate in the Milky Way disk of 6 = 2.65 per century. 
The total number of SNRs in the disk is then given by 



^ = 61 rsN : 



(14) 



for rsN = 1-83 10^ years we reproduce the value given by 
Leahy & Xinji (1989), who estimate iV = 485 ± 60 for 
Galactic SNRs with a limiting surface brightness of E > 
3 X 10^22 Wm^2 Hz^^ sr^^. From Eq. ( [l3| ) we obtain: go — 
1386 ± 171. 



The difference is due to the fact that with the older value of Rq 
10 kpc in Eq. ( px| ) the peak of the distribution appears at rmax 
~ ~ ~ 4 kpc 



a o 



5 10 15 20 

Galactocentric Radius 

Fig. 1. Surface density of SNRs and supernovae, respectively, 
in galaxies as a function of Galactocentric distance R. The 
solid curve, Q{R), gives a fit to the SNR distribution of the 
Milky Way and is given by Q{R) = 87A5R^-^ exp[-0.322i?], 
whereas the dotted curve is a fit to the extragalactic supernova 
distribution by Dragicevich et al. (1999), given by S{R) = 
350 exp[— 0.185i?]. Note the excellent agreement in the expo- 
nential decrease of Q{r) and S{R) for 5 < i? < 20 kpc. 

It should be emphasized that the SNR lifetime, tsn, does 
not enter the analysis of Sect. ^ (see Eq. (pl|)), since it is only 
the CR production rate and hence the SNR rate per unit volume 
that determines the truly diffuse 7-ray emissivity. 

6. Flattening of the radial 7-ray emissivity 
gradient 

In the following we will argue that the CR propagation pic- 
ture needs to be changed for two physical reasons. Firstly, as 
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we have already pointed out, diffusive transport in the disk and 
in the halo is not the same due to a different magnetic field 
geometry and to different wave excitation and damping mech- 
anisms. Thus anisotropic diffusion seems to be rather the rule 
than the exception. Secondly, advection of CRs is the dominant 
transport mode above a certain transition boundary located at 
(Rc, zc), with zq — zc{r) varying with Galactic radius. 

More specifically, we consider SNRs as the primary sources 
for the Galactic CRs. Each remnant results from an energy de- 
position of £^sN — 10'^^ erg, which is converted into relativis- 
tic particles with a certain efficiency i' ^ 0.1 — 0.5 (Drury et 
al. 1989; Berezkho & Volk 2000); the individual value of ly 
is poorly known from theory due to uncertainties in the injec- 
tion process and depends e.g., on the value and orientation of 
the circumstellar magnetic field. SNe exploding inside a su- 
perbubble, i.e. a hot tenuous medium, initially generate low 
Mach number shocks, which are less efficient in accelerating 
particles. After a time of the order of the sound crossing time, 
however, the shock impinges on the much colder and denser 
surrounding shell and becomes progressively stronger thereby 
accelerating particles more efficiently; this should to lowest or- 
der compensate for the initially decreased efficiency in the hot 
ambient medium. We would expect that due to the continuous 
energy input by successive SN explosions also a long lived 
shock would be able to accelerate particles to energies in ex- 
cess of 10^^ eV, albeit adiabatic energy losses would become 
more and more severe with time. Since the diffusion coefficient 
of CRs increases with energy, advective transport of particles 
significantly above 1 GeV will be the dominant mode of trans- 
port only at larger distances from the Galactic plane. Therefore 
these particles will quickly fill an extended halo and not gener- 
ate many 7-rays in the disk via TTp-decay. 

Bykov & Fleishman (1992) have argued that successive ex- 
plosions inside a bubble can generate strong turbulence, which 
should transform a significant amount of the total free energy 
to cosmic rays. However, at the same time the injection rate 
at the shocks may be reduced as a result of shock modifica- 
tion due to previously generated CR particles. With the details 
of the acceleration mechanism in superbubbles being still de- 
batable, we believe that to lowest order there is no difference in 
the overall energy transfer from thermal plasma to CRs, if a SN 
explodes inside a superbubble or just forms a single remnant. 
Thus the energy production rate of CRs should be roughly pro- 
portional to the number of SN explosions, regardless whether 
they occur in the general ISM or inside a superbubble. The nu- 
merical difference in the derived CR energy density (and CR 
pressure) (cf. Eqs. (|l8|) and ([l7|)), however, between treating 
particle acceleration in superbubbles as equal to single rem- 
nants, and disregarding acceleration in superbubbles altogether, 
is small. According to the previous Sect, it amounts to a factor 
^tot/^SNR = 2.65/2.2 w 1.2, and is therefore well below the 
uncertainty in the acceleration efficiency i^. In the following, 
we tend to be conservative and retain a low value of ly = 0.1. 

Using the results from Sect. ^ we can estimate the local (in 
radius) number of SNRs within a circular ring of the Galaxy 
with a width of, say, 2 kpc. Relating the SN rate directly to 



the number of observable SNRs by Eq. (|l^ and writing (cf. 
Eq. (^)) 



Qi{x) 



AN{x) 



(16) 



27ra;Aa; 

we can now define the SNR rate per unit radial confinement 
volume, rj{x), by 

AiV(x) 

ri{x) = 



TSN 

dQijx) 
NzcRl 



(17) 



The local CR energy density, ec in the disk at z w 0, can 
then be estimated to be 



7c 



7c - 1 



Pc{x) = i rj{x) EsN V Tcsc , 



(18) 



with /SV{x) = 2ttxAxzcRq, and Tcsc being the local dif- 
fusive CR storage volume and time, respectively; the factor i 
accounts for one half-space, e.g. z > 0. The CR pressure in the 
disk is denoted by Pq (x) and 7c ~ 4/3 is the specific heat ra- 
tio for the CRs. Working in the framework of an average tran- 
sition boundary, zq, which divides space below and above in 
purely diffusion and advection regions, respectively, the prob- 
ability for particles at z > zc to return and hence contribute to 
the CR pressure in the disk becomes exponentially small. We 
have therefore used z = zq as the vertical extension for the 
CR storage volume. As we can see from Eq. ([T7|), ri{x), de- 
pends on the local and total number of sources and the global 
SN rate, but not on the SNR life time. The vertical extent of the 
diffusion-advection boundary zc is determined by the balance 
of diffusive and advective CR flux at this vertical distance, i.e. 
to lowest order by 



D 

Zn ^ — 
^ V 



D 



uo + Vaq ' 
and the local CR escape time is 

zlix) D 



D {m{x)^VAa{x)f 
Thus the radial CR pressure variation is given by 

7c - 1 ri{x)vEsNZc 



(19) 



(20) 



27c {uo[x) + Vm{x))' 
7c - 1 9Qi{x)vEsN 



(21) 



27c N RI{uo{x) + Vm{x))' 

It is easy to show that Eq. ( pi] ) follows directly from the kinetic 
equation for CRs, as can be seen from Eq. (|47|). 

The diffuse 7-ray intensity resulting from Tr'^-decay pho- 
tons is 



(X I nc n-adl ^ (nc) (nn) I , 



(22) 



with (nc) and (rtn) being the line-of-sight averaged CR 
and hydrogen number density, respectively, and / = 
■y/z^ + (r — Rq)^ is the line of sight. From Eq. (|T]) we have 

Qi{x) 



nc (X Pc{x) oc 



Uo + Vaq 



(23) 
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disregarding the changes in the CR distribution function due 
to energy-dependent diffusion, nn is an observationally fixed 
quantity, for which we have to use its line of sight averaged 
value (see Eq. (p2[)) if we are interested in quantitative calcu- 
lations of the 7-ray intensity; here we are primarily concerned 
with the radial variations of I^. 



In order to obtain uq 
we should assume that 

(uo + Vao) oc Qi(x) , 



const for a weak 7-ray gradient, 



(24) 



if ~ const; if V — Vqz, we have zq ^ \/D/Vo and there- 
fore should require 

{uo + Vao) 'xQlix). (25) 
Combination of Eqs. ( pO[ ) and ( p4[ ) leads to 
D 



oc 



(26) 



All other things being equal, a locally enhanced CR production 
rate should therefore lead to a reduced CR escape time and 
diffuse 7-ray intensity. 

7. Radially Dependent Cosmic Ray Driven 
Outflows 

As we have deduced in the last Sect., the CR pressure Pc{x) 
in the disk is a radially dependent quantity, and therefore we 
expect the outflow velocity, u{x), and the mass loss rate, m{x), 
to be also radially dependent. In an earlier paper (Breitschwerdt 
et al. 1991) we have shown that such a behaviour already exists 
as a consequence of the radial dependence of the gravitational 
potential. The net result was a monotonic decrease (increase) 
of terminal velocity (mass loss rate) with increasing Galactic 
radius for a radially constant mass density, po- Now we have 
superposed the radial variation of the CR source density Qi (x) 
and investigate in the following how this changes the outflow. 

However, Eq. ( ^H ) is an implicit equation, since Pc{x) 
depends on uq{x) + Vaq{x), which in turn depends on the 
energy density available in CRs, i.e. ec{x) and thus Pc{x) 
itself, to drive the outflow. To that end we have performed 
self-consistent galactic wind calculations of the fully nonlin- 
ear equations, in which for a given gravitational potential of 
the Milky Way, and a relativistic CR gas (7c = 4/3), to- 
gether with a spatially averaged mass density g{x) ^ qq = 
1.67 X 10~^^gcm~'^, an average thermal pressure Pg{x) ^ 
PgO = 2.76 X 10^^'^ dyne cm^^, an averaged halo magnetic 
field B{x) ^ Bq ~ 1/iG, and a small average level of wave 
amplitude 5B{x) / B{x) ~ SBq/Bq = 0.1, the advection- 
diffusion boundary zc{x), Pc{x) and uo{x) + VAo(a^) are cal- 
culated self-consistently, using qq — 1492, the value derived 
from the numerical integration of Eq. (|l2|) (see Sect. ||). The 
form of the potential (including, disk, bulge and dark mat- 
ter halo) and the opening of the flux tube due to geometri- 
cal divergence are the same as used by Breitschwerdt et al. 
(1991). We have used a radial grid for our calculations in 
the Galactic radius interval from 0.5 — 12 kpc, with spacings 
of 0.5 kpc. The net radial dependence of the outflow veloc- 
ity we obtained from our numerical calculations is shown in 



0.4 X Ve,, 

0.5 X (Uio-H Vaio ) 

uo + Vao 
[km/s] 

0.1 X Q(r) 




Fig. 2. Escape speed, Vesc — V— 2<I>o, of the Galaxy (dotted 
curve), outflow velocity, uio + waio. at a vertical distance of 
z = 10 kpc (long dashed curve), and outflow velocity, uq+waoj 
at reference level, zq (solid curve), respectively, as seen in a 
fixed frame of reference, and the supernova distribution func- 
tion (short dashed curve), Q(r), (see Eqs. (|l0|) and ([l3|), as a 
function of galactocentric radius r. Curves have been scaled 
down by the factors indicated on the vertical axis to fit into 
one diagram. For boundary conditions of self-consistent galac- 
tic wind calculations see text. 

Variation of Referenoe Level with Galoctio Rodius 



•^ 2 4 6 8 10 12 

Galactic Radius [kpc] 

Fig. 3. Dependence of the diffusion-advection boundary, zc (r), 
on Galactocentric radius r. The location of zq at any radius is 
defined by the balance of diffusive and advective CR fluxes. 
The curve shown here has been obtained from self-consistent 
galactic wind calculations (see text). 
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Fig. ||, and its r- and z-dependence in Fig. 0. The outflow ve- 
locity mq + Vao in Fig. ^ (solid line) for a fixed Galactic frame 
of reference at the radially varying transition boundary ZQ{r), 
shows a similar convex shape as the source distribution Qi{x) 
(short dashed curve), although much less pronounced. This is 
due to the fact, that the decrease of the gravitational poten- 
tial, as can be seen from the escape speed, v^sc — \/— 2"I>o, 
with $0 = + $d + for bulge, disk and halo contri- 
butions, respectively, (dotted curve in Fig. keeps the out- 
flow velocity from dropping off too rapidly with increasing 
Galactocentric distance. Therefore, even a uniform SNR dis- 
tribution would enforce a radially varying outflow velocity, as 
has been shown in Fig. 7 of Breitschwerdt et al. (1991). The 
radial variation of the outflow velocity at a distance of z = 10 
kpc from the disk, uio + Vaio. is shown by the long dashed 
curve, which both shows the transition from the flow velocity at 
reference level zq to the escape speed, and also a roughly l/r^- 
dependence, which we will use in our analytical treatment later 
(cf. Eq. (H)). 

It can be seen from comparison of zc{r) and UQ{r) + 
Vp^i^{r), as obtained from the fully nonlinear calculations, that 
the simple ansatz of Eq. ( [T^ is indeed fulfilled. The functional 
dependence of zc(r) with radius is straightforward to under- 
stand. Close to the Galactic centre, the gravitational pull is 
strongest, as can be seen from Vcsc in Fig- @; since we chose 
all other quantities being the same across the disk (constant 
density, thermal pressure, magnetic field strength), the outflow 
velocity, and hence mass loss rate, are smallest here. Eq. ( p^ 
then tells us, that zq (r) must be largest. At the outer parts of 
the Galaxy, the gravitational field becomes weaker, but also the 
source distribution, and hence the CR pressure, decrease expo- 
nentially, and so the outflow velocity (cf. Fig. ^ also decreases, 
and zc{r) must increase again (see Fig. It is noteworthy that 
the maximum of UQ{r) + VAo(r) ^nd the minimum of zc(r) at 
r ^ 6 kpc do not coincide with the maximum of Qi{r) and 
Pcir), respectively, at r ~ 3 kpc. This must be a consequence 
of the interplay between the gravitational field and the source 
distribution in the fully nonlinear equations (cf. Breitschwerdt 
etal. 1991). 

Finally we mention that the choice of constant boundary 
conditions across the Galactic disk is a conservative one. In 
reality we should also enhance the thermal temperature and 
pressure in regions of higher supernova activity. The net effect 
would be a more pronounced peak in outflow velocity and a 
deeper minimum in zc(r), respectively, and therefore an even 
"better" quantitative proportionality between UQ{r) + VAo(r) 
and Q\{x) according to Eq. (|4|). 

8. Models of advective and anisotropic diffusive 
CR transport 

In the following we discuss several advection models (some in- 
cluding also anisotropic spatial diffusion), in which we exam- 
ine in detail the ideas presented in Sect. ^. We show that for a 
given radial CR source distribution function oc Qi(r), 

we can find a function f2{r) oc V{r), for which the effect 
of the radially varying Galactic wind velocity leads to an al- 
most uniform CR distribution in the disk. We refer to relations 



/2(r) ~ F{fi{r)) as "compensation equations". It is thus pos- 
sible to explain the observational data by pure CR propagation 
effects. 

To fix ideas, we start out with a functional relation for the 
Galactic CR source distribution of the form 

Qiz,r,E)^K{E)hir)S{z), (27) 

neglecting the thickness of the Galactic disk; here S{z) denotes 
the Dirac delta function. Then for a given source distribution 
function, /i (r), we try to find a wind structure of the form 

V{r,z)^V{z)hir), (28) 

which makes the CR distribution almost uniform. In the follow- 
ing we describe models of increasing complexity with respect 
to geometry and outflow velocity. 

8. 1. One-dimensional model with a constant 
galactic wind velocity 

We start from the simplest one-dimensional CR transport 
model, which is an extreme case of anisotropic diffusion, since 
only propagation in the z-direction is allowed. In this case the 
galactocentric radius r is simply a parameter in the model. 

The one-dimensional diffusion-advection equation for CR 
nucleons can be written as 

where for simplicity we take the advection velocity to be con- 
stant and directed away from the Galactic plane 

F(r,z) = ±Fo/2(r)e(±(z-zc)). (30) 

Here 8(z — zq) is the Heavyside step function, and the sign 
"+" corresponds to the velocity above, and the sign "— " to the 
one below the Galactic plane. 

The boundary conditions for CRs are determined either by 
free escape into intergalactic space, if the density of electro- 
magnetic fluctuations generated by the CR flux decreases fast 
enough (see Dogiel et al. 1994), or by CR advection in a galac- 
tic wind to infinity (Breitschwerdt et al. 1991; Ptuskin et al. 
1997), if the level of fluctuations is high enough. Which of 
these cases is relevant for the Galaxy, is the subject of a separate 
investigation. Fortunately, CR spectra and densities in the disk 
are independent of the boundary conditions /ar away from the 
Galactic plane, if there is an outer region of advective transport. 
In the case we discuss here, it is assumed that the CR prop- 
agation region can be formally divided into a diffusion halo 
wrapping around the Galactic plane and an adjacent advection 
region, reaching out to intergalactic space. These two regions 
are separated by a boundary surface, zc, at which both the CR 
density and the flux have to be continuous. 

The sources are concentrated in the disk and are supposed 
to emit a power-law spectrum of particles 

QiE,r,z)^ZdKh{r)E-'"S{z), (31) 

where Zd is the thickness of the Galactic disk and the coefficient 
K determines the CR production by SNRs in the disk. 
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The location of the transition boundary for a constant dif- 
fusion coefficient D and advection velocity Vq is determined in 
this approximation by 



zc 



D 



(32) 



Then for z < zq the transport equation reads 



D- 



dz^ 

and for z > zq 



ZdK fi{r)E-^'> 6{z), 



^{ncVof2ir)) + ^Vof2{r)S{z-zc)ho~l)nc^0. (34) 

Since adiabatic energy losses do not change the particle spec- 
tral index, n cx E '"' . The solution of Eq. (|3|) is 

nd = Ci + C2Z (35) 

and the CR density is constant in the advection region, 

nc - Co . (36) 

The constants can be determined from the boundary conditions 

at z = 



D^^z^KMr)E-'- 
clz 



(37) 



and at z = Zc by continuity of the diffusive and advective CR 
fluxes 



Since at z — zq also = nc, we obtain 

^ 3zdKE-^o hjr) 
K)(7o + 2)/2(r)' 

and therefore 

rid = ZdKfi{r)E^ 



z- Zc 



(38) 



(39) 



(40) 



.■K)/2(r)(7o + 2) D 

We see that at z = 

z,KMr)E-'o{^o + 5) 

hV, (70 + 2)' ^^'^ 

and the "compensation" equation must have the form 

/iW = /2W, (42) 

for the CR particle density to be uniform, in other words, the 
radial dependence of the outflow velocity must be the same as 
that of the source distribution. 

From Eq. ( p^ we can also estimate the total pressure, Pc, 
of CRs in the disk, which is 



00 

Pc ^ j nd{z = 0, E)E ■ dE 

Eo 

^ z^KhE-^^"-^^ (70 + 5) 
ilo~2)f2Vo (70 + 2) 



(43) 



(44) 



The physicals units of K are [K] = erg''''' ^/s, and 
equals the CR power in the Galactic disk, hence 



KE, 



'(70-2) 



(neglecting constants of order unity) 



N 



(45) 



On the other hand the normalized source distribution (cf. 
Eq. dl^)) is equivalent to 



(33) Q^{x) AN fiAVd 



-flQ 



A^d 



AAd 



(46) 



where and Vd are the area and volume of the Galactic disk, 
respectively, and identifying Vq with uq + Vaq, Eq. be- 
comes 



KE-'^-'^-'^ (70 + 5) 
(70 - 2)f2Vo (70 + 2) 
i(^) 



vE, 



SN- 



(47) 



N{ua + Vao) ' 

which corresponds to Eq. (|2l|), thus verifying the validity of the 
intuitively derived relation between CR pressure, source distri- 
bution and Galactic SN rate in a Galactic wind flow. 

8.2. One-dimensional model with z-dependence of 
the wind 

Based on the conclusions of the previous Sections, we can gen- 
eralize the one-dimensional solution obtained by Bloemen et 
al. (1993), taking into account radial variations of the sources 
and the wind velocity. Let us suppose that the advection veloc- 
ity has the form 



V{z,r)^3Vozf2{r); 

the solution of Eq. is then given by 

n{z = 0, r) ~ A—^^k)=E"' , 
DVof2{r) 



(48) 



(49) 



where A is a constant; note that Vo is a typical time constant 
(e-folding time) of the flow. 

We find that a vanishing radial dependence of the CR dis- 
tribution, n = const, can be obtained if 



f!{r)^f2ir). 



(50) 



In these simple models the size of the halo can be derived 
only locally, because global or large-scale characteristics, like 
e.g. the CR distribution along the Galactic plane, are indepen- 
dent of the CR propagation itself. 

These simple analytical solutions of the one-dimensional 
transport equation show that the main effect, which formally 
leads to the "compensation", is a curved transition boundary 
between diffusion and advection regions, i.e. 



zc{r) 



D 



(51) 



note that such a behaviour is indeed apparent from Figs. ^ 
and H The main point here is that the compensation effect 
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is a consequence of the r-dependence of the CR life-time, 
zl/D ^ D/V'^[r). 

In the framework of the pure one-dimensional model it 
is unimportant how far from the Galactic plane the boundary 
zc{r) is, but as we shall see, this figure will be essential for the 
three-dimensional case. 



8.3. Three-dimensional axisymmetric model with 
radial velocity dependence 

To demonstrate the effect of radial changes of the boundary 
zc{r) between regions of diffusive and advective propagation 
of CRs, we investigate the diffusion-advection equation with 
a given radial dependence of the wind velocity and the source 
density in a more realistic cylindrical geometry (axial symme- 
try, i.e. d/dtf = 0) with the velocity varying as 



r 

and the source distribution as 



Q{i;z) = = Qofiir 



(52) 



(53) 



Then the equation for the CR distribution function we have to 
solve reads 



Vo dn „ 9^ 1 d f dn 

— "5~ = Dz-^n + Kr-— r— 
r az oz'^ r or \ or 



(54) 



where we have stressed the anisotropic diffusion by two dif- 
ferent constant diffusion coefficients in the radial and perpen- 
dicular direction, and D^, respectively. Note that Vb has the 
same dimensions as the diffusion coefficients, and that n{r, z) 
is the particle distribution function in phase space. The bound- 
ary conditions at z = are implicitly included in Eq. (|4|). Here 
we treat a special case of radial variation of the boundary zq- 



z^ir) ~ r—, 



(55) 



corresponding to the velocity dependence given in Eq. (p2[). 
Comparison with Figs. |^ and || shows that such a behaviour 
may indeed represent the outer Galaxy, where the outflow ve- 
locity decreases with increasing radius and the distance of the 
boundary zq from the Galactic disk increases accordingly. 

Eq. (|^ can be solved analytically if we restrict ourselves 
to self-similar solutions, in which the distribution function does 
not vary with r and z independently. The price we have to pay 
is that there is no unique solution that covers both < r < 00 
and < z < 00. Instead a similarity solution is found for 
z ^ and r > 0, and one for z > 0, which have to match in 
overlapping regions. We start out with the latter one, applying 
a transformation of independent variables in the form 



(56) 




and introducing the function 

<^(^) = z^^n{z,r) . 



Equation ( p4| ) in terms of and then reads: 

g2j, a J, 



+(64 + 6500- 

where the coefficients are given by 

Vo 



0, 



(58) 



h = 1, 



62 - 



65 = 



^(^- 1) 



(59) 



This equation can be cast into a self-similar form for two sets 
of the parameters a and fj,: a — 2, n — —1, and a — 1, n — 0. 

We can derive an analytical solution of this equation in the 
form (see Bakhareva & Smirnova 1980) 

0(0 = X(l +^^)'^exp(AarctanO , (60) 

where if is a constant determined by the CR source power. 

The unknown coefficients /i, C and A are determined from 
the system of algebraic equations which one can obtain from 
Eq. ( pSf ) by equating the coefficients for the different powers of 
^ to zero. The resulting system of algebraic equations is given 
by 



Abi + 64 = 

2C{l + bi) + A{A + b2) + b5 = 

2C{2A + 62) - A{2 - 63) + 64 = 

2C(2C - l + fe3) + 65 = 0. 



(61) 
(62) 
(63) 
(64) 



We note that this system for the unknowns A, C and /i is 
overdetermined; therefore we have to check for consistency. 
From the system of equations (|6T1)-(|64|) it can be deduced 

that 



A 

C 
Ml 

M2 



2 ' 
0, 
-1 . 



(65) 

(66) 

(67) 
(68) 



It can be verified that these solutions are all consistent. 
However, /i = implies A = and C = 0, and hence does not 
yield any physical solution. Taking [i — —\ gives C — —\. 
Therefore the solution is given by 

(/) ^ K{1 + ^2)-i/2 exp(AarctanO , 



(69) 



where the constant K is determined from the conservation of 
particle flux at ^ = 0. 

The limit z ^ 00 corresponds to ^ ^ so that 

lim0(O = Mm K(l + f^)-^/^ exp(AeiTcteinf) 



1^0 



(70) 



(57) 



because at large distances from the disk (z 00), the particle 
transport is completely determined by advective transport, and 
the number of particles ejected by the sources is conserved. 
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The full solution of Eq. (Q) reads 

Or, 

n(r, z) = — 



formal integration near ^ = gives the boundary condition at 

1 = 0, 



Vo 



exp ■ 



1/2 



arctan 



iD^r 

Z 



(71) 



(72) 



d$(0) 



^0 



=$(0) = - 



Qo 



(81) 



From the solution Eq. (|72|) we can obtain the particle density 
behaviour at large r and z, which are: 

(i) for z ^ 00: 



The function $ is obtained by solving the homogeneous 
equation 



+ l)^ + [(l-2^)e--^]§ 
The solution of Eq. (K2h is found to be 



= . (82) 



1 



n (X — exp 



(73) 



(ii) and for r 



1 



00: 



C2 



(83) 



14) ^ 

n (X — exp I 



(74) 



2F1 ( 7], 77, 1 + /X + 77; i (1 ± 



where 
1 



(84) 



From Eq. ( ^4| ) it is clear that we are unable to describe the 
particle distribution near 2; = with the self-similar variable ^ ^ 2 \ ~'~ ^^ D k 
(Eq. (|6|)) and the function from Eq. (^7|), since it is impossible 
to satisfy the boundary conditions at ^ = cx). Our goal is to study the dependence of the particle density 

In order to study the radial behaviour of the distribution near z = 0, in other words we have to determine the value of 
function near z = a different self-similar ansatz is used, A^- 



which can be extended down to the spatial region near the 
sources. The similarity variable has the form 



From the asymptotic expansion of hypergeometrical func- 
tions at I 00, we find that the solution Eq. ( p3| ) depends on 
^ as 



Z 



and the transformed distribution function reads 

$(^) = r^'^n{z,r) . 



(75) 



(76) 



lim $(0 cx , 



(85) 



This also follows directly from Eq. (|77[), which for large ^ has 
the form 



Then Eq. (54) becomes 

d2$ 



(1 - 2Ai)^ 



$ = 0. 



(86) 



1 



de 



(l-2A^)^ 



(77) 



and if 

a + /^-2 = -l, (78) 
the RHS of Eq. (ITTI) can also be cast into self-similar form 



^S{z) 



Qo 



It is easy to verify that <!> behaves like a power law in ^, i.e. 
$ cx (87) 
and the characteristic equation 

A(A-l) + A(l-2/i) + ^2 ^0 (88) 
has just one degenerate solution: 

Ai,2-M- (89) 



-AO- 



(79) 



We have already shown that at large z the distribution func- 
tion n has the form n oc 1/z (cf. Eq. ([73|)). On the other hand 



The boundary condition at ^ = immediately follows from it follows from Eq. (^ for ^ ^ 00 
integration of Eq. ( [ttI ) for small ^, in which case the equation 
becomes 



n cx r^$(0 cx r'^^'' = 



(90) 



V d^ y/D^, 



and therefore /i = —1, thus matching the two solutions ob- 

(80) 

tained for different regions in r — z-space. It is clear that by 
the choice of the similarity variable (Eq. (|7^), it is impossible 



We integrate Eq. ( |80| ) along a box, encircling the Galactic plane to describe the distribution function near the Galactic centre 
with i — 0, from —e to +e, and subsequently let e ^ 0. The (r — > 0), but everywhere else in the disk. 
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Since close to the sources, ^ 0, $ ~ const., we derive 
from Eq. @ 

n(r) cx - , (91) 
r 

whereas the source distribution is required to have a radial ex- 
ponent a = 2 as follows from Eq. (j78|), and therefore 

Qir) cc 1 . (92) 

Thus radially dependent advection flattens the distribution 
of CRs in the disk as compared to the source distribution 

/2(r)(xr-\ /i(r)cxr-2 (93) 

and the "compensation" condition, which in this case means a 
partial compensation since n (x 1/r rather than n ~ const., is 
realized in this case if 

/2 = /i • (94) 



Nucleon Distribution Function 




10 

r [kpc] 



Fig. 4. The nucleon distribution function n(r, z) is plotted in ar- 
bitrary units as a function of Galactocentric radius r and height 
z above the plane. The radial and vertical diffusion coefficients 
are k^. — 3 x lO'^^cm'^/s and = 3 x lO^^cm^/s, respec- 
tively. 

In Fig. ^ we have plotted the resulting nucleon distribution 
function, which shows a 1/r-dependence (cf. Eq. ([T^), thus 
flattening out the steeper radially dependent source distribu- 
tion ({X by a radially varying transition boundary. It 
should be emphasized that such a behaviour is a natural con- 
sequence of the underlying radially varying source distribution 
itself. All previous models with zc{r) — const, have ignored 
the growing influence of advective transport by an increase in 
CR sources and hence local CR pressure. Of course, it should 
be kept in mind that zq is also energy dependent, since from 
observations of the CR spectrum it is inferred that D oc E^, 
where the value of S depends on z-variations of the Galactic 
wind velocity. 



Nucleon distribution function 




10 

r [Itpc] 



Fig. 5. Same as Fig. ^ but the radial and vertical diffusion coef- 
ficients are now different, viz. Kr = 3 x 10^^ cm^ /a and = 
3 X lO^^cm^/s, respectively. This corresponds to anisotropic 
diffusion occurring preferentially in the disk, where the mag- 
netic field is oriented mainly in radial direction. 



Nucleon distribution function 




r [kpc] 



Fig. 6. Same as Fig. ||, with Kj. = 3 x 10^® cm^/s and = 
3 X lO^^cm^/s, respectively. This corresponds to anisotropic 
diffusion occuring preferentially in the halo, where the mag- 
netic field is oriented mainly in the vertical direction. 



Moreover, according to Fig. ^ a stronger radial than verti- 
cal mixing of CRs due to anisotropic diffusion leads to a fur- 
ther flattening of the nucleon distribution function and hence 
to a weaker diffuse 7-ray gradient. We expect this to happen in 
the Galactic disk, where the geometry of the large scale mag- 
netic field is mainly parallel to the disk. In contrast, in the halo a 
substantial vertical B-component might exist (like e.g. in NGC 
4631), resulting in a larger value of as compared to k,-. 
As can be seen from Fig. || this effect is competing to some 
extent with the also vertically directed advection velocity, and 
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therefore the distribution function is similar to the one shown 
in Fig. I 

In summary we conclude that the value of the halo exten- 
sion calculated from the gradient data in the framework of the 
isotropic diffusion model may indeed be an artifact. In the next 
Sect, we outline a general solution. 

8.4. General solution for axisymmetric convective 
and anisotropic diffusive CR transport 

We now want to work out a general solution of the stationary 
two-dimensional CR transport equation for nucleons, without 
any restriction of the radial behavior of the source distribution 
Q{r, z). In the following we use iV(r, z, E) for the nucleon dis- 
tribution function in order to avoid confusion with the number 
n for the enumeration of the poles in our solutions. 



d 



V (0a;, E)VN ~ u{x)N) ~ — f -Vu{x) E N 



1. 



(95) 



where the diffusion tensor, D_, in cylindrical coordinates with 
axial symmetry is given by 



D = 



7.7. 



K, 

■ E^ 



(96) 



We are not so much interested in energy dependent than spa- 
tially anisotropic diffusion and advection, and so, for conve- 
nience, we set 6 = 0. For the nucleon component other than 



adiabatic losses are negligible (dE/dt 
axial symmetry reads 

r dr 



dz 



Vo djNE) _ 
~ 72 dE ~ 
with 



Q^2 

Qi{r) 6{z) E-^^ 



-- 0); thus Eq. 

,Vo d{Nz) 

J.2 



(97) 



(98) 



and where we have chosen a spatial variation for the velocity 
field as 



u{x) = 









SVof2{r)z 



(99) 



and a power law spectrum for injection by the disk sources. 
The choice of u{x) is motivated not only by simplicity, but 
also by our numerical calculations (see Figs. (|]) and (Q)). The 
latter figure shows the outflow velocity as observed in a fixed 
Galactic frame of reference, V{r, z) = Uz{r, z) + Vx{r, z). It 
can be seen that V{r, z) increases roughly linearly with z and 
decreases in a power law fashion with r, so that Eq. ( ^9| ) is a 
good first approximation. The boundary conditions for obtain- 
ing these results are given in Sect. ^. In future calculations, we 
will also test other analytical approximations to V{r, z). 

It is clear from Eq. ( |97| ) that the effect of anisotropic diffu- 
sion corresponds to stretching or compressing scales in respec- 
tive directions, since in the new variables, r = r/y^ , and. 



V(r, z) 
[km/s 




z [kpc] 



r [kpc] 

Fig. 7. Galactic wind outflow velocity (z-component), 
V{r, z) = Uz(r, z) + vj^{r, z), as a function of Galactocentric 
radius r and vertical distance z from the Galactic plane, 
according to self-consistent Galactic wind calculations for a 
given SNR distribution (for details see text). 



z = z/vDz, we recover the equation for isotropic diffusion. 
Obviously, this may also change the gradient of the cosmic ray 
density in the disk, but we will show that rather strong modifi- 
cations of the density (including the observed uniform distribu- 
tion of cosmic rays) can be obtained for spatially nonuniform 
advection. 

In order to solve Eq. ( |97| ) analytically, we apply the follow- 
ing transformations 

Nir,z,E) = E-">S{tg), (100) 
Vo 



A = ^ 



B = J — 



7 
f 

Q 




(101) 



(102) 

(103) 
(104) 

f = eB. (105) 
i = ze-<^. (106) 

The differential operators for the variable transformation 

(r, z) —f {g, ^) are given by 

d_ 

dr 
d_ 

dz 

Qj.2 

dz^ 



dg 



(107) 
(108) 



Ad 



d£,dg 



dg^ 



d_ 

dg 



de ■ 

Using the above transformations, Eq. (pTh can be cast into 



, (109) 



(110) 



.2^^''S 



d^S 
dg^ 
Qi(e^) 



2^ 



d'^S 
d£,dg 



(l-3A)e 



dS_ 



(Ill) 
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Next we transform the boundary conditions appropriately: 



N{r,z = ±00) = 
N{r +00, z) = 
dN{r = 0,z) 



dr 



= 



±00) = (112) 
S{g= +00,0 = (113) 
'4) =0. 014, 



as J.. 



In order to proceed with the solution of Eq. (Ill), we look for 
the Green's function G(^, g) of 



- A^G = b{q- d)m- 

with boundary conditions 



G{e,i^±^) 



(115) 



(116) 
(117) 



A standard method for solving PDE's of the above kind is to use 
integral transforms; since we consider a steady state problem 
over all r- and 2:-space (or ^ and q, respectively) we use Fourier 
transforms of the kind 



G{i, g) 



—F{k,^)exp{~ikg) 

-00 



G*(^, g) ex.p{ikg)dg . 



(118) 
(119) 



Inserting the transform Eq. (118) into Eq. (115), and after inte- 
grating by parts we apply the boundary conditions Eqs. ( |116| )- 
( |1 17[ ) and end up with an ODE in ^, 

(1 + e)F"{0 + i'2^k-3A+l)^F'{0 

- (fc2 + A-f) F{0 = S{0 exp{ikg') . (120) 

Let us inspect the associated homogeneous equation first. 
Eq. (120) can be cast into the hypergeometric differential equa- 
tion 



x{x - l)y" + [{a + P+ l)x - 7] y' + aPy = , 

by applying the following transformation 

X = p + iq^ 
=4> dx = iqd^ , 



and 

x{x - 



(121) 



(122) 
(123) 



(124) 



from which we infer p — 1/2 and q = ±1/2. Matching the 2nd 
term of Eq. ( 1 2 1 ) we obtain 



Solving for a and /3 from Eqs. (126) and (127) gives 
1 



{2ik -3A±G) 



P1.2 = -{2ik-3ATG) 



(128) 
(129) 



where C = y^A{9A - 12ik + 47). In Eqs. ([128]) and ( [l29| ) a 
and P are complementary and we can drop the subscripts, using 
"H-" for a and " for f3. 

Thus the solution of the homogeneous equation reads 



F{0 = C12F1 (a,/3,7;-(l±*C)) + 



C2 



(i±iO 



1-7 



(130) 



2F1 ( a-7+l,/3-7 + l,2-7;-(l±iC) 



where 7 is not an integer and < Re{l ± < 1. 
Applying the linear transformation 

2Fi(a,/3,7;a;) = (1 - x)^""-'^ x 

2F1 (7 - a, 7 - /?, 7; x) , 



(131) 



we finally obtain 



FiO = C12F1 (a,/3,7;-(l±ze) ) + 







1-7 




C2 






■^(IT.^)" 



- 7 — a— /3 



2F1 ( l-a,l-/3,2-7;i(l±zO 



- C12F1 ( a,/?,7;-(l±zC) ) + 



C2 



1-7 



. 1 



2F1 ( l-a,l-/3,2-7;-(l±zO 



(132) 



where the integration constants ci and C2 have yet to be deter- 
mined from the bound ary c onditi ons ( see Appendix From 
the Fourier transforms (US) and (y_9) we see that this implies 



F(^ = ±oo) = G; 
furthermore we require 

(i) continuity of F(^) across ^ = 

(ii) = FiO . 



(133) 



a + (3 + 1 = 1-3A + 
7 = l:{a + f5 



2ik, 
f 1) 



(125) 
(126) 



We also note from Eq. ( 132 ) that F(^) is multivalued in the 
complex plane and therefore a branch cut has to be introduced, 
running along the imaginary axis from —i to +i. 



Note that there is no difference between the cases with "+" and 



We now return to the inhomogeneous equation Eq. (120) 



-". The last term in Eq. (121 ) requires 



For ^ 7^ 0, Eq. (1133) is a valid solution; in order to evaluate 



a/3 = -(fc2 +7v4) , 



(127) 



the solution at ^ = we integrate Eq. ( 12C ) in the infinitesimal 
range — £ < ^ < £, with |£| ^ 1 and then take the hmit £ — > 0. 
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Exploiting the properties of the solution F(^) outlined and using the expression (136) for 0) — F{£^ — 0), we 
above, and noting that the integration over an odd function van- obtain 



ishes, we obtain after a little algebra 

F'{Q+)-F'{0-) = exp{ikg'), 

which is a jump condition for the derivative at ^ 
Integrating this result over ^ gives 

F(0+) -F(0") = 0, 



(134) 
= 0. 

(135) 



N{r,z^ 0,E) 
1 



•-\-ioo r 



Qi{r')2Trrdr x 

f)r(V) 



r(i±^)r 1- 



(143) 



as expected. With Eqs. ( |133| ) a nd (|134| ) we can n ow evaluate the 
constants ci and C2 (see Eqs. ( B.12 ) and ( B.13 ), Appendix^. 

Manipula ting t he F- functions with help of the relations in 
Appendix ^ ( |B.9| - [B.11[ ) the solution for ^ = reads 



(- 


1)T'C0S 


(¥) 


+ sin 


V 2 ) 


(- 


l)T'sin 1 




+ cos 


f TTa \ 

V 2 / 



-K 



dK . 



m 



zexp(zM) f^(t)r(f)r(i-^)r(i 



With the help of the Green's function from Appendix |^ (cf. 
Eqs. ( |C.1C| ) and ( |C.11| ) we obtain 



27r2 



1 



xsin(f) 


cos(^) 




^cos(f) 


+ sin(^y 


(-1)7 sin (f) 


+ cos(^) 



(136) 



p+00 

/ Qi{r')2TTr'dr' 

J r 

(-l)^cos(^) +s 



r(f)r 1^ 



8.5. Solution in tiie disk 



Using the solution for F{£_ = 0) we can now, as a special case, 
calculate the full solution in the disk, i.e. N{r, z — 0, E)). The 
Green's function is then given by Eq. (118) 



G(e = 0, Q- d) 



+00 fHi- 

—F{k,0)exp{-tkg). 
27r 



(137) 





t)r(V) 



rmr(i-f) 

ffy^dK 



(144) 



In our solution (136) we have suppressed the explicit depen- 

dence on fc, but of course it is implicitly contained in a(fc), ^j^^^^ integration contours ^ and Cp correspond to r' > 



(3{k), and7(fc) 

Let K = ik, we obtain 



and r' < r, r espe ctively. The only poles to consider for A 3> 1 
are (see Fig. C.l ) for n = 0, 1, 2, ... 



G{0;g,g') = / F{K,0) exp[-K{g- g')]dK a(i^„) = -2n =^ K,, = -2n ± A{-f + 6n) , (145) 



where F{K, 0) = exp{-K g')F{K, 0). Then 



where the sign "— " should be chosen for Kn- The poles a{Kn) 
are inside the contour Ca (see Appendix^. 

Then we get from the theorem of residues (the residues 



G{^; g' y^^ ' e^' 2TTe<^' dg' . (139) of Gamma functions r(n) at n = 0, -1, -2, ... ai-e given by 



i-ir 



) 



Transforming back to the variable r and r' (Eq. (105 )): 

GiO;ry) = -l^FiK,0){-) dK , 
and 



N{r,z = 0,E) 



+ 00 



G(e;r,r') 



hn2^(rr^dr'. 
D,. ^ ' dr' 



(140) 



(141) 



AnD, 



-E 



+ 00 



E 



Qi{r' )2t: r' dr' 

(-i)"r(i^) 



The particle distribution in the disk is then given by 



N{r, z:^Q,E) = 



2'kD, 



-E- 



^ n!r (i^) r (1 - %) 

(_l)7.cos(^) 
(_l)7„sin(^)+(-l)'^ 



(146) 



(7) 



dKn 

da 



-\-ioo 



Qi{r')2Trr'dr' F{K,0) (^^'^ 



-K 



dK (142) 



Here 

/3„ 



-2n -iA- 2^/1(7 + 671) 



(147) 
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1 



7n 

da 



1 - 4n - 3A - 2^^(7 + 6?- 



= 1 



-2n - y/A{^ + 6n) 
3A 

2^^(7 + 6n) ■ 



(148) 
(149) 
(150) 



Although Eq. (147) looks rather awkward, this solution for 
A 3> 1 exhibits simple asymptotic dependencies at different 
radii. If the CR sources occupy a limited volume of the disk, 
bounded by a radius a, Q{r') — Qor'^'^Q{a — r'), then for 
r > a the function N is zero, since the integral J^°^ = 0. 
This means that advection can remove particles from the plane 
z = so fast that none remain in the plane at radii r > a (see 
Appendix^). 

Inside the source region (r ^ a), N{r) can be written as 



+00 



7V(r)cxrV^/ r'^'-'-^^dr, 



(151) 



so that for strong advection A ^ 1 the integral is determined 
by the lower limit independent of the source distribution. Then 
we have 



N{r) cx r 



2-q 



(152) 



and the CR density is almost constant in the disk region if (7 « 
2. We have therefore demonstrated that in the case of strong 
and spatially nonuniform advection, the CR density in the disk 
exhibits only a weak gradient in spite of strong radial variations 
of the source distribution. Such a behaviour may be responsible 
for the mild gradient observed in the diffuse 7-ray data. 

In the diffusion dominated case { A <C 1) the values of the 
integrals are determined by the source distribution. Here we 
should take into account the integral over Cp, which is not zero 
(see Appendix Far away from the disk, i.e. r ^ a, the dis- 
tribution function is determined by the only pole, (3 — \, which 
is inside the contour Cp (see Fig. |C.5| ). Therefore independent 
of the source distribution, we have 



N{r) cx r-i-V'4(7-3) ^ ^- 



(153) 



just what is expected for pure diffusive particle propagation. 

Therefore we conclude from these analyses that a more or 
less uniform CR distribution in the disk can be expected, if 
advection is strong, i.e. A ^ 1. 

The "compensation" equation in this case reads 



(154) 



The interesting conclusion from the solution of the general 
three-dimensional model is that we can reach a quasi-uniform 
CR distribution especially when the boundary zc is closest 
to the Galactic plane in regions where the source density and 
hence advection is strongest (see also Fig. i.e. when A ^ 1. 
This is in sharp contrast to the results of the pure diffusion 
model, because only an extremely extended galactic halo can 
produce a more or less smooth CR distribution and thus a mild 
diffuse 7-ray gradient. 



9. Source contributions to thie 7-ray bacl<ground 

The examples of the last Sect, substantiate our argument that 
transport effects should be responsible for all but eliminating 
the 7-ray gradient in the GeV energy range. This ignores possi- 
ble contributions to the 7-ray flux from the sources themselves. 
In fact, SNRs as CR sources are expected to accumulate the ac- 
celerated CRs before releasing them into the ISM at the end of 
their lifetime. In nuclear collisions with thermal gas atoms or 
(for electrons) with photons inside the remnant, this high den- 
sity of energetic particles will produce a strong 7-ray intensity 
that contributes a "diffuse" background if these sources are un- 
resolved. In fact the ensemble of Galactic SNRs, expected to be 
strongly concentrated in the disk, will constitute such a back- 
ground with typically less than about 10 such sources along a 
line of sight in the disk. They will be unresolved in GeV 7-rays, 
except for a few well-known nearby objects (Berezhko & Volk 
2000); also their collective flux contributes less than 10% to the 
overall 7-ray flux. However, we know that the CR source spec- 
trum is much harder cx £;-(2 ito2.2) jj^^j^ jj^g diffuse Galactic 
CR energy spectrum (cx E~'^ '^), and therefore the source con- 
tribution will equal the truly diffuse 7-ray flux at about 100 
GeV and become even dominant at about lO'^ GeV (= 1 TeV). 

If the sources collectively dominate in the form of an unre- 
solved "diffuse" background at very high 7-ray energies, then 
it should be possible to observe the average CR source distri- 
bution and even the average CR source spectrum (!) along the 
Galactic disk in the TeV range. Present observations at TeV 
energies have not yet been able to detect this "diffuse" back- 
ground. However, the experimental upper limits and the pre- 
dicted "diffuse" flux are so close to each other (Aharonian et 
al. 2001) that an instrumental increase by an order of magni- 
tude should lead to a detection. The next generation of TeV 
instruments like H.E.S.S., CANGAROO, MAGIC or Veritas, 
that has this sensitivity, will come on line soon. The compari- 
son of such a detected TeV - distribution with the radial shape 
of the source distribution from radio studies would constitute a 
stringent empirical test on our theoretical arguments about the 
transport effects on the truly diffuse Galactic CR component. It 
also shows the basic importance of 7-ray surveys over a large 
range in energies and radial distance. 



10. Discussion and Conclusions 

In the detailed analysis of the previous Sections it was shown 
that a physically reasonable CR transport model should include 
the nonlinear dynamical effects of the Galactic distribution of 
the CR sources themselves. Taking into account a simple radial 
dependence of the source power as suggested by radio obser- 
vations of supernova remnants and pulsars, and thus averaging 
over all azimuthal angles, leads already to a radically different 
spatial behaviour of the CR distribution function. In all cases 
we have studied so far, there is a clear tendency to smooth out 
peaks due to enhanced particle injection by the sources in the 
ensuing radial particle distribution. This can be attributed to 
CR advection in an outflow aided or even caused by the CRs 
themselves. In other words, the importance of CR advection is 
proportional to the CR power of the underlying sources. As a 
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consequence, e.g. the CR escape time, which is usually inferred 
from secondary radioactive isotopes, need no longer be a glob- 
ally averaged quantity, but could vary strongly from place to 
place in the Galaxy. Therefore, the usual assumption that pa- 
rameters of the CR flux measured near Earth are representative 
for the global processes of CR propagation in the Galaxy may 
be a misleading oversimplification. 

It is common practice to determine from the intensities of 
different CR components near Earth the average diffusion co- 
efficient in the Galaxy, the velocity of advection, the height of 
the CR halo in the direction perpendicular to the Galactic plane, 
and the CR injection spectrum, just to name the most important 
ones. Based on this hypothesis the nearly uniform radial CR 
distribution, derived from the measurement of diffuse Galactic 
7-rays, can be reproduced only if there exists thorough spatial 
mixing of CRs in the framework of an extended halo (if CR 
diffusion is isotropic). Hence in such a case local and global 
properties of CRs do not differ from each other. However, the 
inferred halo height from chemical composition (zq — 2 — 4 
kpc; see e.g., Bloemen et al. 1991, 1993; Webber et al. 1993; 
Lukasiak et al. 1994) is clearly inconsistent with the value de- 
rived from the interpretation of the 7-ray data {zq = 15 — 20 
kpc; cf. Dogiel & Uryson 1988, Bloemen et al. 1993) within the 
framework of an isotropic diffusion model (see Appendix 
We therefore conclude that the halo size derived from CR nu- 
clear data reflects only a local value near Earth, and the huge 
halo extension derived previously from 7-ray data may be an 
artifact, since it relies on the validity of global values for lo- 
cally obtained CR data. This conclusion is supported by our nu- 
merical galactic wind simulations, which show that the vertical 
distance of the diffusion-advection transition boundary from 
the Galactic plane, is inversely proportional to the CR source 
power and not spatially constant as been previously assumed. 

Radio observations of external galaxies indicate a large- 
scale magnetic field geometry, which is mainly parallel to the 
major axis in the disk, and if a halo field exists, it is parallel 
to the minor axis. Therefore we expect that CR diffusion is in 
general anisotropic, with a radial diffusion coefficient Ki in the 
disk, which is much larger than diffusion in the perpendicular 
direction, D^, and vice versa in the halo. In this case the ini- 
tially inhomogeneous CR distribution, due to a radially varying 
source distribution in the disk, is smeared out, whereas in the 
halo the dominant diffusion component can be superposed 
by a strong advection velocity, which may determine the spatial 
particle distribution. 

It would be desirable to have a high enough spatial reso- 
lution and photon statistics in the future to observe the radial 
distribution of diffuse 7-rays above 100 MeV in nearby edge- 
on galaxies, such as NGC253. However, it seems unlikely that 
both space-borne and ground-based 7-ray observatories will 
satisfy this requirement in the near future. Thus the only direct 
observation of the CR source distribution in the Galaxy will be 
possible with next generation TeV instruments like H.E.S.S. 
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Appendix A: The halo size from iocai 
characteristics of the CR flux 

The flux of CRs observed near Earth (locally) contains in- 
formation about regions from where these high energy parti- 
cle have been produced, as well as about the location of the 
boundary from where these particle escape into the intergalac- 
tic medium and do not return to the galactic disk. This means 
that in principle we can infer the size of the diffusion region 
simply from local CR properties. Let us denote the thickness 
of this region by Zosc- We will demonstrate how to derive this 
quantity from a set of simplified equations, describing densi- 
ties of stable and radioactive secondary nuclei generated by 
primary CRs in the gaseous disk. 

For a one-dimensional diffusion model the system of equa- 
tions for the density of stable. Tig, and radioactive, n,-, nuclei is 
given by 



D 



dz^ 
and 



dz^ 



^NgVa,pZdnld{z) 



(A.l) 



(A.2) 



The distribution function of the nuclei in the disk is not very 
sensitive to the boundary conditions. Therefore we analyze 
here the simplest case, in which at Zosc it is required that 



Us, Tlr = ; 



(A.3) 



here aSp and (Tip are the spallation cross sections of stable and 
radioactive isotopes, respectively. We assume that processes of 
nuclear spallations take place in the gaseous disk with half- 
thickness Zd, which is much less than the halo scale height ^esc- 
Here z is the coordinate perpendicular to the Galactic 
plane; the gas is supposed to be concentrated in the Galactic 
disk with z^ <C ^osc, where the gas density equals A^'g, and rip 
is the density of primary cosmic rays in the disk (z = 0). Then 



AT / 

Us = Ngvasp-rip— (zcsc 
and correspondingly 



z) 



7A\/riT sinh( ^°^'= 
and thus the ratio ^ at z is 



■) 



/ = 



tanh 



(A.4) 



(A.5) 



(A.6) 
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The CR life-time in the region of diffusion propagation is de- 
termined from 



CSC 

IT 



and for r ^ Tcsc we obtain 

rj„ /T. 

/ 



sp 



rp 



■'sp •^csc 



(A.7) 



(A.8) 



We see that this ratio depends on the two unknown parame- 
ters D and z^sc- Therefore they cannot be determined indepen- 
dently. To eliminate one parameter, we need additional obser- 
vational constraints, for instance the average thickness of mat- 
ter traversed by CR nuclei in the Galaxy during their lifetime, 
given by the expression 



2(J 

X = MpNgV Tcsc , 

^csc 



(A.9) 



where v is the velocity of the nuclei and Mp is the proton rest 
mass. 



Combining this with Eq. (A. 8) we obtain 

2 



X ~ 



MpZiiNcrVT ( U-i: 



-'sp 



(A. 10) 

-12- 



which gives Zcsc — 1 kpc, for z^Ng ~ 4 • lO^^cm ^, x 
14 g cm~2 (Ferrando et al. 1991, Ferrando 1993), r = 2.2 • 10^ 
years and / ~ 0.25 (see e.g., Simpson & Garcia-Munoz 1988). 
Here we took <Trp ^ (Jsp for the estimate. 

If we include the convection term into the kinetic equation 
the analysis becomes more complicated, but the result is almost 
the same. If the velocity depends on z as 



Viz)=Woz, 

we obtain for the ratio / (see Bloemen et al. 1993) 
/ 



or, 



(A.ll) 



(A. 12) 



which requires Vq — 10 km s ^ kpc ^. The grammage in this 
case is 



nN„V , 

P « D{E) ' 



(A.13) 



where ^;esc is the coordinate of the surface above the Galactic 
plane from where CRs are sucked away by the Galactic wind 
(i.e. the transition boundary between diffusion and advection). 



IPjE) 
Vo 



(A. 14) 



fsp) 



Combining Eqs. (A. 12) - (A.13) we also obtain for (Crp ^ 

(A.15) 



MpZdNgVT 



f ^csc 

and Zosc is also about 1 kpc. 

This value estimated from the chemical composition of 
CRs near Earth is usually defined as the size of the CR halo. 
On the other hand the halo size can be derived from the global 



7-ray emissivity distribution in the Galactic plane. The value of 
Zosc then determines the size of the CR mixing volume, which 
is estimated to be at least 10 kpc. If the diffusion model would 
describe the CR distribution correctly and satisfactorily in the 
Galaxy, the estimates of Zosc obtained from local and global 
data would be similar. However, the sizes of the CR propaga- 
tion region differ from each other by about one order of magni- 
tude, ruling out pure diffusion as the sole CR transport process. 



Appendix B: Determination of thie integration 
constants cl and c2 

Using the transformation formula 



r(/?)r(7 - a) 

2F1 ^a, 1 — 7 + a, 1-/3 + a; 

r(7)r(a-/3) , 

r(a)r(7-/3)^ 



2F1 (/3,l-7 + /3,l-a + /3;- ) ,(B.l) 



observing that x = i(l±-g^), and the limit (which of course 
will also hold upon interchanging a and (3 due to their symme- 



try, see Eqs. (126) and (h27b) 



1 , (B.2) 



we obtain from the condition ( |133| ) for Eq. ( |132[ ) after some 
algebra 





Atto. 




— a 


lim < 




X 




[ 







2F1 q;,1 — 7 + a,l — /3 + q;;- 

X 



Cl 



r(7)r(/3-a) 
r(/3)r(7 - a) 



C2 e 



~27r{7 — 1) 



1) r(2-7)r(/3-«) 



2 V 



r(/3-7 + l)r(l-a) 

/3 



2F1 /?, 1 - 7 + /3, 1 - a + /3; - 

X 



Cl 



C2 e 



r(2 - 7)r(a - /3) 



r(a-7 + l)r(l-/3) 



r(7)r(«-/3) 
r(a)r(7 

= . (B.3) 



Observing that Re{a) > and Re{b) < 0, Eq. ( |B.3[ ) can only 
be satisfied in general, if 



(B.4) 



D. Breitschwerdt et al.: The gradient of diffuse 7-ray emission in the Galaxy 



19 



Next we use condition ( 1 34 ) to fix ci and C2 . Making use of the Eq. (132) becomes 
relations 



dx 



2F1 (a,/3,7;a;) = 
a (3 



— 2Fi(a + l,/3 + l,7+l;a:) , 
7 



Fi a + l,/?+l,7+l;- = 



1\ _ 40F7 r(7) 



ai3 r(f)r(f) 



2-a,2-/3,3-7; 



4V^r(3-7) 



(l-a)(l-/3)r(W(i^) 



(B.5) 
(B.6) 

(B.7) 



we obtain 



lim F'(f) 

5^0 



exp(ifcg') 



±47r- 



2B 



Cl 



r(7) 



r(f)r(§) 



+ C24T 



^-1) r(2-7) 



r(i^)r(V) 



(B.8) 



Inserting Eq. (B.4) into Eq. ( |B.8| ) we end up after some tedious 
algebra and using the following relations 



r(l + z) = zT{z) 

r(z)r(i-z) 



r(2z) 



sin(7rz) 



r(z)r z 



(B.9) 
(B.IO) 

(B.ll) 



with 



Cl 



and 



C2 



i_B exp(iA;p') 



(-irr(f)r(f)sin(i^) 



2\AF 



r(7)[(-l)^sin(^)+cos(^)^ ' 

(B.12) 



exp(iA:e') 2-("+« cos {^) T{^)T{^) 



= 0) = T 



i ex:p{ikg') 



r(f)r(^) 



^ 2 



^(-ircos(f; 


+ sin(^)^ 


^(_l)7sin(f ) 


+ cos(^) 



(C.2) 



Then the Green's function for ^ = reads 

G{C = 0,g)=J _F(fc,e = 0)expHM 

If K = ik, then the transform for the variable r is 



FiK,0)i^-j dK, 



-K 



(C.3) 



which is determined by the residues of the poles of the F- 
functions 



and 

P{K,n) = 2to + 1 , 

and by the residues of the poles determined from 



(-l)''sin 



tt(3 
"2" 



+ cos 



(t) 



0. 



(C.4) 



(C.5) 



(C.6) 



Since 7 = a + /3 + 1 this equation can be written in the form 











r. / /? 




exp 


ITT \a + — ~ 




— exp 


^n (a--- 






exp 
= 0, 
yielding 

/3 = -a-l + 4n 



(-.-i.f)] 



+ exp 



(C.7) 



(C.8) 



However, these values of f3 also make the numerator vanish, 
i.e. 



r(2-7) [(-l)^sin(f)+cos(^)^ 

(B.13) 



(-ircos(^)+sin(:^)=0 



(C.9) 



Appendix C: Determination of integration 
contours for tlie Green's function f, r') 



The Fourier transform 

^+00 



because all we are doing is to interchange a and /?. Thus (3 = 
—a — 1 + 4n does not give any poles. 

We see from Eq. {Z.3) that the contour of integration should 
be enclosed in the semiplane Ke{K > 0) for f > r' , i.e. the 
contour in Figs. Z.l and |C.2|. In this case the integration 



G(^, g) cxp{ikg)dg , 



over the outer part of the closed contour is zero and the func- 
tion G equals the residues of the sum inside the closed contour 
^^■^^ only. For f < r' the contour is enclosed in the outer part of the 



semiplane Re(if < 0), i.e. the contour C/3 in Figs. Cl and C.2 



of the Green's function of Eq. ( 11 5 ) is shown in Eq. (119) with The locations of the poles in the -fC-plane are determined from 
the constants ci and C2 taken from Eqs. (B. 12 ) and (B. 13 ). For a one-to-one map of the contours onto the a- and /3-planes and 



the special case of the solution in the Galactic plane (^ = 0) subsequently Eqs. {CM and (C.5) are applied 
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The Green's function reads 



k - plane 



X da 

da 

dKp 
X — - 

d(3 



dp \ , (C.IO) c 



and applying the theorem of residues we obtain 

G(^-0,r,r') = 





Rek 



Fig. C.l. The contours of integration Cq and Cp in Eq. ( C.1CI| ) 
for f < r' and f > r' in the fc-plane for 1. To show these 
contours in the K = ik plane the figure should be rotated by an 
angle of 7r/2. The dashed circles represent poles of the second 
kind, the soUd circles poles of the first kind (see text). 



(C.ll) 



(-l)T'sin 



where a, /3, 7 and K are determined by their values at the poles 
of the Gamma functions. The residues of the function T{n) 



have the value 



(-!)'■ 



at n = 

The contours for f > 



0,-1,-2,. 
r' and f < 



in the fc-plane are 



shown for A'3> I in Fig. C. 1 , and for A ^ 1 in Fig. C.2 



To determine the locations of the poles of the F-functions 
a{Kn) and (3{Km) in the X-plane we should map a unique 
array of K values one-to-one onto corresponding arrays of a 
and (3 values. In this way we define the contours of Re(_ft') > 
and Yie{K) < in a and /3 planes. 

We start with r > r' . We find from a (if = 0) (see Fig. C.3) 
that the contour Cq intersects the abscissa at the points: 



and 



iA 



1 



47) 



(C.12) 



Therefore all poles q;„ are inside the contour Cq, . On the other 
hand the area encircled by Ca does not contain poles We 
see that poles mapped onto the fc-plan e are transformed to 
the negative part of the real axis (see Figs. ^ ^ and 
Since there is a branch point in the fc-plane at 




Rek 



(C.13) 



Fig. C.2. The contours of integration Ca and Cp in Eq. ( ClOj ) 
for f < r' and f > r' in the fc-plane for A <C 1. To show these 
contours in the K = ik plane the figure should be rotated by an 
angle of 7r/2. The dashed circles represent poles of the second 
kind, the soUd circles poles of the first kind (see text). 



the contour of integration Cp (r > r') has a cut (see Figs. C.l 
and |C.2| ). Its reflection in the /3-plane is shown in Figs. |C.4| and 
|C.5| . We see that the contour Cp intersects the abscissa axis at 
the points 



and 



■iA 



1 
3 



(C.15) 



3A 
2 



1 



yjAiQA + 47) 



(C.14) 



The value of the latter is negative for A ':$> 1 and positive 
for A < 1 (7 ~ 4). Therefore only if A < 1 the pole 
/3 = 1 is inside the contour Cp (see Fig. 2.5 ). Correspondingly, 
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a - plane 




-6 -4 -2 
^ @ @ @- 



Im a 



Re a 



Fig. C.3. The contour of integration Ca for r < r' in the em- 
plane. The sohd circles represent poles inside the contour Cq. 



A Im (3 



(3- plane 




Fig. C.5. The contour of integration Cp for r > r' in the j3- 
plane for A 



ther e is o nly one pole inside the contour Cp in the A; -plane (cf. 
Fig. C.2). In t he op posite case, A':^ \, all poles are outside this 
area (see Fig. C.4) and then the integral over Cp equals zero. 
This means that the poles jSm mapped onto the A; -plane are o f 
the second kind (their phases are larger than 2tt, see Fig. C. 1 ). 



The reason, is that in the case of strong advection and the ve- 
locity increasing towards the center as r~^, advection removes 
cosmic rays from the central part of the disk so fast that these 
particles do not reach its outer parts. Therefore the contribu- 
tion of the integral over Cp (the central part of the disk) is zero 
and the Green's function is determined by the integral over Cq 
only. 



A Imp 




Fig. C.4. The contour of integration Cp and the poles /3„ 
f > r' in the /3-plane for A^ 1. 



for 



We conclude this Sect, by deriving the transformations be- 
tween K-, a- and /3-planes, and by deriving the Green's func- 
tions both for strong and weak advection. From the condition 
of one-to-one correspondence we obtain for the poles in the 
if-plane 



Kn ^ ~2n - y^Ai-f + 6n) . (C.16) 

• K,-,-, = 2m + 1 + 



Since 

aiK„) - -2n 
P{K^) = 2m + 1 , 

for n,m = 0, 1, 2, ... we have 

/3„ = -2n-3A 



\/v4(7-3(2m + l)).(C.17) 

(C.18) 
(C.19) 



6n) 



7n = ^ (1 - 4n - 3^1 - 2v/^(7 + 6r 



Kn = ~2n ~ y/A{j + 6n) 
dK„ _ 3A 



da 

7m 



1 



2^yl(7 + 6n) 
2m+l-3A + 2V^(7- 3(2m+ 1)) 

i(^2(2m + l) -3.4 + 1 + 
2V^(7-3(2m+l))) 



Km = 2to + 1 + v/A(7 - 3(2m + 1)) 
dKra ^ ^_ 3A 



dp 2^A(7-3(2m+l)) 
and for A > 1 (cf. Eq. (|C.10|)): 

G(C = 0,r,r') = 



(C.20) 
(C.21) 
(C.22) 
(C.23) 

(C.24) 

(C.25) 
(C.26) 
(C.27) 
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r °° (-i)"r 

e(r'-r)^ 



t nlT (i^) r (1 - 



(C.28) 



_(_l)7„sin(^)+(-l)" 
For A <C 1 the Green's function reads 

G(C = 0,r,r') = 

r <x, (-i)"r 

S e(/-.)i: 




da 



l-lir 



da 



X 



";ton!r(i^)r(i-^) 

(_l)7.sin(^)+(-l)"| 

f sin (^) I /^x-^^ dif, 
|(_l)7^+cos(^) J Vr' 



d/3 



where 

a/3 

7/3 



d/3 



1 - 3A + 2 V^(7 - 3) 
i(2-3A+l + 2V^(7-3)) 



= 1 



1 + v/^(7 - 3) 
3A 



2v'A(7 - 3) 



(C.29) 

(C.30) 
(C.31) 

(C.32) 

(C.33) 

(C.34) 
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